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We consider the quantum XY model and study the effects of interacting perturbations on the time 
evolution of the von Neumann and Renyi entropies of spin blocks after global quenches. We show that 
the entropies are sensitive to perturbations that break hidden symmetries behind the integrability of 
the model. At times much larger than the characteristic time of the well-known linear increase of the 
entropies, we identify a time window characterized by a novel linear growth followed by saturation. 
The typical time of the phenomenon is inversely proportional to the perturbation strength and the 
behavior is trigger off by the extinction of an infinite number of local conservation laws following 
a non-abelian algebra. The universality of the crossover is revealed by a semi-classical picture that 
captures the leading behavior of the entropies. We check our theoretical predictions against iTEBD 
simulations. The good agreement between theory and numerics substantiates the method developed 
in [Bertini and Fagotti, J. Stat. Mech. (2015) P07012] for investigating a pre-relaxation limit in 
weakly interacting models. 


I. INTRODUCTION 

A global Hamiltonian parameter is suddenly changed 
and the state, originally in equilibrium, evolves under 
a different Hamiltonian. This simple protocol of non¬ 
equilibrium time evolution is generally called a “global 
quench”. Since the first studies on global quenches, many 
exciting theoretical aspects have been uncovered (emer¬ 
gence of statistical descriptions, universality of dynam¬ 
ics, etc.)^^^^ and observed in numerical simulations^"^^^^ 
and real experiments^^^"^^. One of the first theoretical 
insights was that the entanglement entropies (von Neu¬ 
mann and Renyi entropies) of subsystems follow a uni¬ 
versal behavior^. In the framework of conformal field 
theory (CFT), Ref. [1] showed that the entropies of an 
interval grow linearly in time until a particular time pro¬ 
portional to the subsystem length i. Then, they saturate 
to stationary values. It was soon realized that the lin¬ 
ear behavior is almost independent of the system details: 
despite CFTs describing only a class of critical models, a 
similar time evolution appears in the (space-time) scal¬ 
ing limit i ^ t ^ 1 also in noncritical spin chains. The 
semiclassical interpretation proposed in Ref. [1] is indeed 
expected to hold rather generally^^^^®. Furthermore, the 
extent of the time window of linearity is inversely pro¬ 
portional to the (Lieb-Robinson^®) maximal velocity vm 
at which information propagates, so the analysis of the 
entropies represents a simple way to extract vm- 

After Ref. [1] there has been an increasing interest in 
the non-equilibrium time evolution of the entanglement 
entropies. A large part of works concerned the dynamics 
induced by short-range Hamiltonians^'*’'*®’^^^®®, but re¬ 
cently some attention moved also to global quenches with 
long-range interactions®^’®^. In addition, the entangle¬ 
ment entropy has been at the center of intensive research 
in other non-equilibrium problems, like local and geomet¬ 
ric quenches®®”®®. It is also worth mentioning that there 
are a few proposals, based on local quenches, to measure 
the entanglement entropies in real experiments®^’®®. 

In spite of the numerous studies on global quenches. 


not much attention was paid to the relaxation process of 
the entropies: the lack of a simple picture explaining the 
approach to the stationary values has probably been a 
major deterrent to the research on that topic. Indeed, 
to the best of our knowledge, the long time behavior of 
the von Neumann entropy after a global quench has been 
analyzed analytically only in the quantum XY model^®, 
where, depending on the system details, it was found 
either a simple power-law relaxation ~ H? jt or a more 
elaborate decay ^ log(</^)/t®. 

The goal of this paper is to identify some universal 
aspects that emerge in a long time limit. We choose a 
backward point of view, asking firstly what are the fea¬ 
tures that dramatically affect the values of the entropies 
at infinite time after the quench. 

It is widely believed that at late times after global 
quenches in generic models local degrees of freedom can 
be described by a Gibbs ensemble (GE) with an effective 
temperature depending on the initial state®®”®”^. In in- 
tegrable models, instead, infinite information about the 
initial state is retained and the stationary properties 
are described by a generalized Gibbs ensemble (GGE) 
constrained by infinitely many (quasi-)local conservation 
laws^’®®. Among the integrable models it is also con¬ 
venient to distinguish the ones with an infinite num¬ 
ber of local conservation laws that satisfy a non-abelian 
algebra®®, which will be refereed to as non-abelian inte¬ 
grable models^^. Also in the latter case a GGE is ex¬ 
pected to describe the stationary properties of local ob¬ 
servables, but infinitely more information than usual is 
retained about the initial state at late times. Glearly, the 
entanglement entropies are very sensitive to constraints 
on the dynamics, and indeed relaxation to a GE, a GGE, 
or a “non-abelian” GGE result in different asymptotic 
values. 

Our strategy is to perturb the Hamiltonian in such 
a way that the limit of infinite time does not commute 
with the limit of zero perturbation. In other words, we 
consider perturbations that break hidden symmetries be¬ 
hind the integrability of the model. It is then reasonable 
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to expect that, as long as the time is sufficiently large 
but much smaller than some characteristic time deter¬ 
mined by the perturbation, the entropies settle close to 
the stationary values associated with the unperturbed 
integrable model (prethermalization plateau): the GGE 
of the unperturbed model, possibly including also non¬ 
commuting charges, is a good local approximation of the 
state. Whenever the unperturbed set of local charges 
is abelian and the perturbation breaks integrability, at 
larger times the expectation values move from the pre¬ 
thermalization plateau and reach the values predicted 
by the GE. When instead the perturbation breaks non- 
abelian integrability into integrability, a subset of charges 
is destroyed and the remaining ones are deformed; as a 
consequence, the earliest plateau (pre-relaxation plateau) 
is generally followed by a less constrained GGE (although 
persistent oscillatory behavior could remain). If the per¬ 
turbation breaks also integrability, in even longer times 
all charges but the Hamiltonian become extinct, inex¬ 
orably leading to a GE. Therefore, we recognize the 
strength of a perturbation that breaks integrability or 
non-abelian integrability as a relevant parameter to the 
late time behavior of local observables and entanglement 
entropies. 

Since infinitesimal perturbations are sufficient to trig¬ 
ger off a crossover between macroscopically different en¬ 
sembles and the entropy per unit length is robust under 
a large class of transformations, we can expect most of 
the details of the system and of the perturbations to be 
qualitatively irrelevant: the time evolution of the entan¬ 
glement entropies in the time window of the crossover 
could display some form of universality. 

a. The model. We consider the Hamiltonian of the 
XYZ model 

= E > (1) 

e 

where (t“ act like Pauli matrices on site £ and like the 
identity elsewhere. We have not explicitly written an 
overall multiplicative constant with the dimensions of an 
energy, which can be viewed as implicitly attached to 
the time. The XYZ spin-^ chain is an interacting in¬ 
tegrable model that can be solved by algebraic Bethe 
ansatz®®. For A = 0 it is reduced to the quantum XY 
model and can be mapped to a noninteracting chain of 
spinless fermions® ‘ . The XY model has an infinite non- 
abelian set of local conservation laws^^, so, in the limit 
A —)■ 0, (A 2 ) describes a non-abelian integrable model. 
We also consider two perturbations that break integra¬ 
bility (for generic 7 and A): 

= and S^=^-Y,a!. ( 2 ) 

t i 

Thus, in its most general form, the Hamiltonian reads as 
H = Ha + A{US^'^ + hS^). (3) 


We multiplied the terms that (are supposed to) break 
non-abelian integrability by A, so A is expected to char¬ 
acterize the timescale at which the hidden symmetry 
breaking becomes manifest. On the other hand, for h or 
U different from zero, relaxation to a thermal ensemble 
generally occurs on much longer timescales that depend 
on At/ and Ah in a different way (the results of Ref. [ 68 ] 
suggest timescales ^ {AU)~‘^,{Ah)~‘^). 

In Refs [22] [69] some effects of hidden symmetry break¬ 
ing have been identified in the time evolution of local 
observables, provided that the initial state breaks one- 
site shift invariance, like the ground state of Ha in the 
antiferromagnetic phase. A posteriori, this explains why 
the original investigations^® on the time evolution of the 
entropy after quantum quenches in the XY model did 
not show any peculiar behavior beyond the scaling limit 
of Ref. [1]. We therefore consider initial states in which 
the symmetry with respect to translations by one site is 
(spontaneously) broken. 

b. Results. Let us denote by 

p, = Tr^P(t))(vI/(t)l] (4) 

the reduced density matrix of a subsystem A consisting 
of £ adjacent spins {A is the complement of A) and by 

(„) _ logTrp“ 

" 1-a 


the corresponding Renyi entropies (the von Neumann en¬ 
tropy being Si = limQ,_,,i+ The main result of this 

paper is that, up to corrections that approach zero in the 
limit of small perturbation, the entropies per unit length 
asymptotically behave as follows 
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( 6 ) 


Here vm is the maximal velocity at which information 
propagates and vm is another velocity emerging at inter¬ 
mediate times; crg°'^ is the entropy per unit length in the 
earliest plateau while a[°‘^ is the value reached at larger 
times. The last undefined inequality signifies that the fi¬ 
nal regime is not expected to extend over infinite times. 
In particular, can be still different from the value 
approached at infinite times after the quench. 

The first time window in ( 6 ) corresponds to a well- 
known limit® and will not be discussed here. On the 
other hand, times t ^ A“® have not yet been explored, 
so we will focus on the limit A <C 1 with At ~ 0(A®). 
We will give a semi-classical interpretation to the scaling 
law ( 6 ) and demonstrate the formula both numerically 
and, in some cases, analytically, by deriving an asymp¬ 
totic analytic expression for Si/£ {ef. (59)). 

c. Organization of the paper. Section H is a concise 
review of what is meant by non-abelian integrability and 
what are the expected effects on quench dynamics. In 
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Section III we summarize the results of Refs [22] [69], 
which are focussed on the dynamics in an intermedi¬ 
ate time window called pre-relaxation limit. Section IV 
provides a physical interpretation of ( 6 ), based on very 
general arguments. The subsequent sections are instead 
devoted to the explicit solution of the dynamics in the 
model with Hamiltonian (Al). Specifically, Section V 
describes the initial state. Section VI is devoted to the 
analytical and numerical analysis of the pre-relaxation 
limit after quenches from the Neel state. In Section VII 
dynamics characterized by persistent oscillations are con¬ 
sidered. We summarize our results in Section VIII. Sec¬ 
ondary aspects and some details of the iTEBD simula¬ 
tions are reported in the supplemental material [70]. 


II. NON-ABELIAN INTEGRABILITY 

The notion of integrability for a quantum model like 
(A2) is deeply connected with its conservation laws. Us¬ 
ing algebraic Bethe ansatz'^, one can indeed construct 
an infinite family of localconservation laws Qj in invo¬ 
lution with one another 

[Qj,Ha]=0, [g.,Q,]= 0 . (7) 

In principle, there could be other independent 
(quasi-)local charges^^^^®, however, in the absence of 
manifest symmetries (like SU{2) for 7 = 0 and A = 1), 
the set of local conservation laws is generally assumed to 
be abelian (7). 

As a matter of fact, there are exceptions in which hid¬ 
den symmetries generate an infinite number of additional 
local conservation laws, the simplest example being the 
quantum XY model, with Hamiltonian Hq^^. The XY 
model can be mapped to free fermions and the dispersion 
relation satisfies Sk = Efc-i-ir • This in turn produces an in¬ 
finite number of local conservation laws Qj that break 
one-site shift invariance and do not commute with one 
another 


[Qj,Ho] = 0, i[Qi,Qj] = fijkQk ■ (8) 

We stress that, relaxing the requirement of locality, the 
existence of noncommuting conservation laws is not ex¬ 
ceptional. For example, any reflection symmetric spin 
model that can be mapped to noninteracting fermions 
has a dispersion relation with the symmetry Sk = £-k- 
This generates an infinite number of noncommuting 
charges, which however are nonlocal. On the other hand, 
the non-abelian set of charges of the XY model is local 
and infinitely larger than the abelian set found in the 
absence of the symmetry Sk = Sk+m justifying the ap¬ 
pellation of “non-abelian integrable model”. 

The immediate consequence of ( 8 ) is that in the ther¬ 
modynamic limit not all local charges are needed to clas¬ 
sify the states. Any maximal abelian subset of {Qj} can 
be used to that aim, and indeed the classification is not 
unique but depends on the subset. 


A. Local charges in noninteracting spin chains 


Spin-i models that can be mapped to noninteracting 
fermionic chains can be investigated with specific tech¬ 
niques that, to some extent, go beyond the standard 
Bethe ansatz solution of the model. Here we briefly re¬ 
mind the reader of the free-fermion formalism to con¬ 
struct the local conservation laws. We focus on one-site 
shift invariant models that are mapped to noninteract¬ 
ing fermions by the Jordan-Wigner transformation (here 
written in terms of Majorana fermions {a^,a„} = 2d in) 

a2i-i = a2i = . (9) 

j<t j<e 

This means that, apart from boundary terms (which are 
irrelevant, if neglected both in the mapping from spins 
to fermions and in the reverse mapping to the spins) the 
noninteracting Hamiltonian H can be written as 

1 

H ~ — 'y ] aiHijUj , (10) 

ij = l 

where "H is a block circulant matrix with 2 -by -2 blocks 


U 


( 1 ) 

in 


(T~L2l-l,2n-l T~L2l,2n-l 
\ n2e- l,2n n2i ,2n 

= i (11) 


with £, n S {1,..., L}. The sum is over inequivalent mo¬ 
menta of the form , with n integer (in the thermody¬ 
namic limit the quantization rule is not important). The 
2-by-2 matrix is usually called symbol. If s is a 

divisor of L we can also define the s-site representation 
of the symbol, which is obtained by gathering 
together blocks of s sites (reducing thus the effective size 
of the chain L —>■ j at the price of increasing the local 
space): 

= ^ ( 12 ) 

k\e'^^L/s — i 


Here £,n £ {!,...,L/s} and is a (2s)-by-(2s) 

matrix with the following properties: 

•H('’)t(/c) =-H(*)(fc) n(-)t(^k) = -H^-^\-k). (13) 

To be explicit, ni+ 2 s{i-i),j+ 2 sin-i) = (fi-in)hj- 

A conservation law Q that can be mapped to nonin¬ 
teracting fermions by (A12) and that is translation in¬ 
variant by s sites has a symbol Q^^\k) that commutes 
with 'H^‘^\k). The conservation law is local if Q^’^'^ik) has 
a finite number of nonzero Fourier coefficients. Remark¬ 
ably, if the spectrum of ii^^\k) is degenerate, one can 
generally find more than 2 s linearly independent sym¬ 
bols Q^^^k) commuting with ?i^®)(fc). Any time that 
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this happens independently of k the Hamiltonian has an 
infinite non-abelian set of local conservation laws, in¬ 
finitely larger than the set that arises when is 

non-degenerate. We refer the reader to Ref. [22] for a 
more extensive discussion. 

To make the notation simpler, since we are going to 
consider the time evolution of two-site shift invariant 
states, in the rest of the paper the superscript identify¬ 
ing the representation will be omitted in the two-site 
representation (s = 2) of the symbols: Q^^^k) —>• Q{k). 


B. Stationary properties after a global quench 

The additional conservation laws can play a crucial 
role in non-equilibrium time evolution. The initial state 
somehow selects the maximal abelian subset of {Qj} that 
will be relevant^^. The mechanism is essentially the fol¬ 
lowing. Noncommuting conservation laws generate de¬ 
generacy in the spectrum. There are “privileged” bases 
in which the Hamiltonian is diagonal and in each degener¬ 
ate energy level a single state has a nonzero overlap with 
the initial state: the dynamics are essentially the same 
as in the non-degenerate case, provided that the Hilbert 
space is properly reduced. From this point of view, we 
can identify the relevant conservation laws as those that 
are diagonal in a “privileged” basis. 

Using the standard arguments behind the emergence of 
stationary behavior®di^ it is reasonable to expect that at 
late times local degrees of freedom will relax to stationary 
values that can be described by a GGE 

PGGE = . (14) 

Zj 

Here are Lagrange multipliers fixed by the integrals of 
motion and Qj are relevant charges, which belong to an 
abelian subset of {Qj}- 

U,Qj]=0 Qj=Y,A,tQi. (15) 

i 

The real rectangular infinite matrix A depends on the 
initial state, which indeed has a twofold effect, determin¬ 
ing both the Lagrange multipliers \j and the charges Qj. 
Alternatively, one could express the ensemble in terms of 
the whole set of linearly independent local conservation 
laws 


PGGE = (16) 

In this way one removes the subtle dependence of the 
operators on the initial state with the drawback of having 
a representation in terms of noncommuting operators. 

So far, non-abelian integrability has simply resulted in 
a (complicated) reduction of the Hilbert space. There are 
however nontrivial effects due to the noncommutativity 
of local charges. We focus on the interesting situation 


in which a small perturbation breaks the hidden symme¬ 
tries behind the additional conservation laws. Roughly 
speaking, two phenomena compete with each other: the 
perturbation tends to approximately preserve a subset 
of the conservation laws, which is however different from 
the abelian set selected (in the sense of (15)) by the ini¬ 
tial state in the absence of perturbations. This gives 
rise to non-trivial time evolution that is disclosed in long 
timescales. 

This phenomenon can be partially understood in the 
framework of stationary perturbation theory in quan¬ 
tum mechanics^® as follows. Let us choose a “privilaged 
basis” of eigenstates of the unperturbed Hamil¬ 

tonian, where n identifies the energy level and i the 
possible degenerate states. The perturbation introduces 
0(A) corrections to energies and eigenfunctions. In non¬ 
degenerate energy levels, up to phases, the eigenfunc¬ 
tions remain “close” to their original values also at times 
t ^ 0(A“^). For degenerate states the situation is in¬ 
stead different: the perturbation splits each degenerate 
energy level (like in the Stark/Zeeman effect) in such a 
way that the energies are approximately equal to the un¬ 
perturbed ones {Enj ^ En'’ + AWnj) but the eigenfunc¬ 
tions mix within the subspace. At times t ^ 0(A~^), 
this results in 0(A°) discrepancies between and 

the corresponding time evolving state: 






|»“(*)) ~ E 


1 ^,( 0 ) 1 ( 17 ) 


3d' 

Here we indicated with Ylj 

the eigenfunctions that diagonalize the full Hamiltonian. 

This analogy fails however in explaining the im¬ 
portance of locality and infiniteness of non-commuting 
charges in the case of a quantum many-hody system. 


III. A PRE-RELAXATION LIMIT 

If the quantum XY model is perturbed so as to break 
non-abelian integrability, e.g. with a term proportional 
to , the expectation values of local observables experi¬ 
ence a crossover^^ in a time window that scales as A“^. 
In generic models, the appearance of plateaux long be¬ 
fore the typical relaxation times is commonly referred 
to as . In Refs [22] [69] 

the crossover at intermediate times was instead called 
pre-relaxation limit, in order to emphasize that it is not 
driven by the breaking of integrability, but rather by the 
breaking of the hidden symmetries behind non-abelian 
integrability; as a result, plateaux at intermediate times 
can be observed also in integrable models. 

For noninteracting perturbations one can rigorously 
show^^ that, in the limit A ^ I with At ~ 0(A°), local 
degrees of freedom are described by a time-dependent 
GGE of the form (16), with time-dependent Lagrange 
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multipliers X£{At): 

p— Ylt Af (At)Qf 

^^- +0(A). (18) 

In Ref. [69] this picture was generalized to interacting 
perturbations and, under a few assumptions, the same 
result has been obtained. The crucial difference from 
the noninteracting case is related to the effective descrip¬ 
tion of the dynamics, which is now generated by a time- 
dependent (noninteracting) Hamiltonian. 


conservation laws and then work out the various terms 
in a mean-field fashion. There are indeed also purely 
nonlocal terms in V, which however are not expected to 
contribute at these time scales®®. The remaining terms 
are polynomials of the local charges, so (21) simply means 
that, in each product of charges, all but one charge, in 
turn, contribute as simple c-numbers. 

We would like to place some emphasis upon the effi¬ 
cacy of the mean-field mapping: quench dynamics in an 
interacting model are reduced to noninteracting (though 
complicated) dynamics! 


A. Mean-field solution 


B. Generalities 


The effective theory proposed in Ref. [69] is like a 
mean-field approximation that takes into account the 
presence of conservation laws breaking one-site shift in¬ 
variance. However, it is not a real approximation, be¬ 
ing instead an emergent description valid in the pre¬ 
relaxation limit. 

Specifically, it was shown that in the limit of small 
A and finite At, the time evolution under (Al) of local 
observables, starting from a two-site shift invariant state 
with cluster decomposition properties^ can be written as 
follows 


( 0 ) (t) ^ ^ (19) 

where the range of the operator must be small in com¬ 
parison with t and 1/A. In (19) (C’)ggEo 
expectation value of the observable O in the mixed 
state p'^^(At), time evolving under a quasi-local time- 
dependent noninteracting Hamiltonian H^^{At) and 
equivalent, at the initial time, to the stationary state 
that emerges for A = 0: 

iO)ZEo {T)=Tt[p^^{T)0] 

zdTP^^iT) = [H^^iT),p^^iT)] ( 20 ) 

pMF(O) = lim lim Tr^[e*'f^'>* l^-o) (^-oj . 

|A|—).oo t—foo 

We introduced the notation T = At to help one identify 
the slowly evolving quantities. The mean-field Hamilto¬ 
nian can be formally written as follows*^ 

/ At/ 

^MF(y) ^ ^ p^{T)Q„ p,{T) = {j^) (T) , (21) 

n \°b;j/GGEo 


where the sum is over linearly independent local conser¬ 
vation laws of the unperturbed Hamiltonian and V is the 
time-averaged perturbation in the interacting picture 


V = lim - 

t—>-(X> t 


n 5A 


^-iHoT 


A=0 


( 22 ) 


The compact notation of (21) means that we take only 
the part of the time average V that depends on the local 


We collect here some useful observations. 

Ensemble representation: 

Since everything is written in terms of the local 
conservation laws of Hq, the mean-field description 
is consistent with (18). Indeed the time evolution 
of a local charge (of the unperturbed model) under 
(21) commutes with the unperturbed Hamiltonian 
(c/. (8)) and hence the mean-field dynamics are re¬ 
duced to time evolving Lagrange multipliers (18). 

Paradox of relaxation in closed systems: 

In the limit 1 ^ ^ <C t the entanglement entropy 
per unit length of a subsystem is generally equiva¬ 
lent to the (thermodynamic) entropy density of the 
GGE describing the late time dynamics*®. Since 
mean-field time evolution is unitary, one could ex¬ 
pect the entropy of the (quasi-)stationary state to 
be conserved. In fact, this is the same paradox aris¬ 
ing in the non-equilibrium time evolution of pure 
states: despite the system being closed, at long 
times reduced density matrices can be described by 
more entangled mixed states [e.g. Gibbs ensembles 
in generic models). Rather counterintuitively, this 
means that if any local observable relaxes in the 
pre-relaxation limit, the mixed state that describes 
the stationary properties has the same form as (18), 
but with a different partition function Z. 

One-site shift invariance: 

It is important to note that the pre-relaxation limit 
for (Al) can be nontrivial only if the initial state 
breaks one-site shift invariance (1-ssi). This can 
be seen as follows. The unperturbed GGE corre¬ 
sponding to a 1-ssi initial state is 1-ssi, being also 
the Hamiltonian (Al) 1-ssi. Analogously, the time 
averaged perturbation is 1-ssi. On the other hand, 
the conservation laws that break 1-ssi are odd un¬ 
der translation by one site®® Qi —>■ —Qi- From 
this it follows that a generic term of the polyno¬ 
mial part of V is generally written as Q1Q2 ■ • • Qn, 
where Qi and Qj take the same sign under 
a shift by one site. By symmetry, the expectation 
value of any odd operator is zero, therefore in (21) 
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only the terms proportional to 1-ssi charges give 
a non zero contribution: the mean-field Hamilto¬ 
nian is one-site shift invariant. We now remind 
the reader that the 1-ssi charges of the quantum XY 
model commute with one another^^. Consequently 
[H^^{T),p^^{T)] = 0 and = p“P(0). In 

other words, the expectation values do not move 
from the values corresponding to the first plateau. 

On the number of non-commuting charges: 

One could wonder whether a finite dimensional 
non-abelian subalgebra of local conservation laws 
could be sufficient to trigger a crossover in the pre¬ 
relaxation limit of the entanglement entropy. In 
that case mean-field time evolution would result in 
oscillations in a finite number of Lagrange multipli¬ 
ers of the time-dependent GGE (16) (the ones cor¬ 
responding to the non-abelian subset of charges). 
This is like a quantum problem with a finite number 
of degrees of freedom, so time evolution is expected 
to be quasi-periodic, with recurrence times depen¬ 
dent on the subalgebra dimension and not signifi¬ 
cantly affected by the subsystem length. The os¬ 
cillations should not modify the extensive part of 
the entanglement entropies, which are in turn ex¬ 
pected not to move from the first plateau (see also 
IV B). As a result, any nontrivial crossover in the 
entanglement entropies can be seen as an effect of 
the unperturbed model having an infinite number 
of non-commuting local charges. 

Problem setting: 

For Hamiltonians of the form (Al), in which the 
interaction terms consist of four fermion operators, 
the (mean-field) dynamics of local observables are 
reduced to the solution of an infinite system of 
quadratic first order differential equations, with an 
extra linear term originated by the noninteracting 
perturbation S^. The system of differential equa¬ 
tions can be derived as follows. First of all one 
has to compute the time average of the pertur¬ 
bation. Then, the purely nonlocal contributions 
are removed and the mean-field differential equa¬ 
tions for the time evolution of the local integrals 
of motion of the unperturbed model are set up. 
This determines the mean-field Hamiltonian in a 
self-consistent way. Being p^^{T) gaussian, the 
expectation value of any local observable can be 
obtained from the two-point fermionic correlation 
function (T^j = Sij — {uiOj)) using the Wick’s the¬ 
orem. Indeed one can show that T is written in 
terms of the integrals of motion (of the unperturbed 
model); therefore, it can be extracted from the so¬ 
lution of the system of differential equations above. 
We point out that, in order to deal with the infinite 
system, the introduction of a cut-off is unavoidable; 
this can be done for example confining the system 
in a fictitious finite box. If the cut-off has been 
chosen in a sensible way the expectation values of 


local observables will be stable under its variation. 

We refer the reader to Ref. [69] for the explicit 
derivation of the differential equations and for a 
preliminary analysis of their solutions. In this pa¬ 
per we will only exhibit the equations in a specific 
example (Section VI). 

Assumptions and checks: 

For interacting perturbations, (19) was derived un¬ 
der two assumptions: the first is similar to the 
hypothesis behind perturbation theory in quantum 
mechanics, in particular the non-diagonal elements 
of the perturbation must approach zero as the cor¬ 
responding unperturbed energy differences vanish. 
The second is that the purely nonlocal part of 
the time averaged perturbation is not relevant for 
the time evolution of local observables in the pre¬ 
relaxation limit (Ref. [69] established that this ap¬ 
proximation is self-consistent). 

We will give grounds for (19) by checking the mean- 
held predictions against iTFBD simulations. This 
will be the most delicate point of this work, being 
difficult to reliably simulate long times. 

IV. PHYSICAL INTERPRETATION 

Let us indicate with rjj the typical range (generally 
independent of time) of H^^{At) (21). 

Ref. [69] pointed out that, in the pre-relaxation limit, 
the system can experience either local relaxation or per¬ 
sistent oscillatory behavior. 

A. Local relaxation 

In the hypothesis of local relaxation, at late (rescaled) 
times At, H^^ becomes independent of time. Goncomi- 
tantly, the expectation values of operators with larger 
ranges start relaxing to stationary values over a timescale 
At that is expected to be proportional to their range i. 
As a matter of fact, this is analogue of what happens 
in the space-time scaling limit I <C t considered in 
Ref. [1]. There, the qualitative behavior of the entropy 
has been understood through a semiclassical reasoning 
that we briefly repeat here. 

Being the energy extensively higher than the ground 
state one, the initial state acts like a source of pairs of 
quasiparticle excitations traveling in opposite directions. 
By translation invariance they are originated everywhere. 
Because the initial state is low entangled, one can as¬ 
sume quasiparticles originated from different points to 
be incoherent. There is also another assumption, namely 
distinct pairs of quasiparticles originated from the same 
point are not entangled. In noninteracting quenches with 
one-site shift invariance one can convince oneself of this 
property by simply writing the initial state in a basis 
that diagonalizes the Hamiltonian. More generally, the 
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FIG. 1. Simplified space-time picture illustrating the scaling 
behavior (6) of the entanglement entropy between an interval 
A and the rest of the system B. The non-abelian integrability 
of the unperturbed model is manifested by the fact that at the 
initial time two pairs of quasiparticles (black and red arrows), 
correlated with one another, with almost identical velocity 
and traveling in both directions are emitted from each point 
(which represents two sites). The light orange region (biggest 
background triangle) shows the subsystems (symmetric about 
the centre) that can be approximately described by the sta¬ 
tionary state emerging in the absence of perturbations (for 
which the velocities of the black and red quasiparticles would 
have been exactly equal). In that regime only the correlations 
between the pairs of particles traveling in opposite directions 
are important (see e.g. the subsystem, horizontal thick line, 
at earlier times). However, at larger times the contributions 
from the cases (highlighted by red circles at the top of the 
diagram) in which a single particle (of the four) crosses the 
subsystem are no more negligible. Such contributions increase 
linearly in time up to the time at which no pair of quasipar¬ 
ticles originated at the same point can cross the subsystem 
together. After that a new equilibrium is reached and the 
subsystems can be approximately described by a stationary 
state that in some cases is the GGE of the full Hamiltonian 
(the green region, namely the smallest background triangle). 
The pre-relaxation limit describes the crossover between the 
light orange region and the green one. 


incoherence of pairs does not seem to be heuristically 
justifiable, but it is in fact appropriate to a large class of 
quenches. Thus, taking that assumption for granted, the 
qualitative behavior of the entanglement entropy can be 
understood in terms of the motion of the quasiparticle ex¬ 
citations. Since the entropy measures the entanglement 
between a subsystem and the rest, it should be propor¬ 
tional to the number of entangled quasiparticles that, at 
a given time t, are both inside and outside the subsys¬ 
tem. This simplified picture explains the linear growth 
followed by saturation observed in the space-time scaling 
limit 1 <C ^ ~ t, at least in the integrable case. 


We borrow this physical interpretation and propose 
some modifications that allow us to understand the be¬ 
havior claimed in (6). Since now translation invariance 
is realized by a two-site shift, each point of the effec¬ 
tive description corresponds to two sites and, in turn, 
the number of quasi-particles originated from the same 
point is doubled. Non-abelian integrability (of the un¬ 
perturbed model) is manifested by the presence of en¬ 
tanglement between pairs of pairs of quasiparticles. This 
is a consequence of the freedom in the choice of the basis 
diagonalizing the Hamiltonian and identifying different 
abelian subsets of charges. The change of basis is indeed 
equivalent to mixing pairs of quasiparticles. 

As depicted in Fig. 1 and explained in its caption, this 
results in linear behavior followed by saturation also in 
the pre-relaxation limit. There are however some differ¬ 
ences with respect to the space-time scaling limit: 

1. Being A ^ 1, at times t ~ the initial state 
can be replaced by the stationary state that lo¬ 
cally emerges when time evolution is generated by 
the unperturbed Hamiltonian Hq, i.e. p'^^(O) in 
(A9) (in Fig. 1, quasiparticles that contribute to 
the crossover in the entropy pass through the light 
orange region). 

2. After the substitution of Point 1, the relevant evo¬ 
lution occurs in a rescaled time At (in Fig. 1, the 
angle between the black and red arrows is 0(A)). 

3. The extent of the linear regime is not related to the 
maximal velocity at which information propagates 
but a new velocity emerges (in Fig. 1, the angle 
between red and black arrows per unit of A is not 
in relation to the direction of the arrows). 

B. Oscillations 

While for noninteracting perturbations like relax¬ 
ation is the norm, it is very common to find persis¬ 
tent oscillatory behavior when the perturbations consist 
of interacting terms. If we insist on the interpretation 
sketched in Figure 1, we realize that oscillations should 
not affect the behavior of the extensive part of the en¬ 
tropy. Indeed we can imagine to deform the quasiparti¬ 
cle trajectories including some oscillatory behavior. The 
amplitude of the oscillations would depend on the quench 
parameters but not on the subsystem length. As a conse¬ 
quence, the effects of oscillations would be concentrated 
on the boundaries of the subsystem, modifying just the 
leading correction 0(£°) to the entanglement entropies. 
This conjecture will be corroborated by a numerical anal¬ 
ysis of the mean-field solution of the problem reported in 
Section VH. 

This is the last section of the first part of the paper, de¬ 
signed to give the reader the opportunity to gain a phys¬ 
ical intuition about the problem under investigation. In 
the following sections the specific non-equilibrium time 















evolution under (Al) will be considered and the mean- 
field solutions worked out in order to derive an analytic 
expression for the time evolution of the entanglement en¬ 
tropies. 


V. THE INITIAL STATE 


Mean-field dynamics are trivial if the initial state is 
invariant under translations by one site (c/. One-site 
shift invariance in IIIB), so we must consider states that 
break such symmetry. We note however that the pre¬ 
quench Hamiltonian can be still one-site shift invariant: 
the symmetry turns out to be spontaneously broken in 
antiferromagnetic or dimerized phases. For the sake of 
simplicity, we choose initial states with factorized form 




cos^ Iti) +e*‘^sin^ |it) 


(23) 


where |t) = It) and |t) = — It). It is not difficult 
to show that, for any value of /3 and cj), l^ko) is a Slater 
determinant for the Jordan-Wigner fermions (A12). De¬ 
spite this property being unnecessary for (19)(A9), it will 
be convenient when, in Section VIB 2, we will introduce 
an intermediate mean-field description. 

There are two relevant states in the class (A4): 


/? = 0 V tt: Neel state (the ground state of iLoo); 


• (/?,()') = (f,7r): product state of singlets (the 
ground state of the Majumdar-Gosh Hamiltonian, 
i.e. (Al) with 7 = 0 , /i = 0, A = 1 and U = I). 

The correlation matrix 


Ffc = 5tn — ('f'olafOnl'I'o) (24) 

of (A4) is a 4-by-4 block circulant matrix. As any block- 
circulant matrix, it is completely characterized by its 
symbol, i.e. the Fourier transform of a block-row {cf. 
(12)). This is given by 

r(A:, 0 ) = cos P — sin j3 {cos -\-sin 4>u'^. 

( 25 ) 

In the pre-relaxation limit the part of the dynamics with 
a typical relaxation time much smaller than 1/A has al¬ 
ready reached (quasi-)stationarity. As a result, the cor¬ 
relation matrix can be replaced by the one corresponding 
to the stationary state for A = 0 . This can be done by 
projecting (25) onto the symbols of the local conservation 
laws of the quantum XY model Hq (A2), which results 
in 


r“^(fc,o) = 

'ycospQsik) - sin/3cos(/)(cos^ |Qi(fc) - qsin^ ^Qsik)) 
cos^ I -I- 7^ sin^ I 

— sin/3sin(/Q4(fc). (26) 


Here Qj{k) identify the eight families of local conserva- 


tion laws of Hq: 



Qi{k) 

= Sk [(T'®e* 2 ‘^ ]( 8 )[cr^e 


Q2{k) 

= cos(fc/ 2 )£fc IG [(T^( 


Qm 

= sin(fc) I G I 


Qi{k) 

= sin(fc/ 2 ) [cr^e*5'""] 

01 

Qsik) 

= £fc [a^e*S‘^^] G [a" 

gie(fc)cr^j 

Qs{k) 

= cos(/c/ 2 ) [cr^e*^'^ ] 


Qiik) 

= sin(fc) G 


Qsik) 

= sin(fc/ 2 )efc [<7^ 


where tan = 7 

tan 1 and Sk = (cos^ 

1 - 1 - 7 ^ sin 


(27) 


IS 

the dispersion relation of the two species of quasiparticles 
that diagonalize the unperturbed model in the represen¬ 
tation in which the local space consists of two sites*®. The 
symbols of a generic local charge can indeed be written 
as cos{nk)Qj{k), with integer n. We point out that for 
j < 4 the associated conservation laws are one-site shift 
invariant (they correspond to the two families of local 
charges identified in Ref. [89]). 

Replacing r(fc, 0) by r’^^(A:, 0) reflects the third equa¬ 
tion of (A9). 


The mean-field description then establishes that the 
noninteracting structure is preserved at 0(A°) up to 
times t ~ 0 (A“^) and the symbol of the correlation ma¬ 
trix remains a linear combination of the symbols of Hq’s 
local charges. 


VI. QUENCHES FROM THE NEEL STATE 


Let us now consider in detail the case /3 = 0 V tt and 
[7 = 0 , namely the non-equilibrium time evolution of 
the Neel state under the Hamiltonian of the XYZ model 
with an integrability breaking perturbation proportional 
to 

The mean-field equations of motion have been worked 
out in Ref. [69] following the problem setting described 
in IIIB. We do not repeat here the (manageable but te¬ 
dious) calculations and simply report the results. 

For states that are reflection symmetric about a site, 
it turns out that the symbol of the correlation matrix of 
p^^{T) (A9) can be written as follows 


r’^^(fc,T) 


E 

ie{2,7,8} 


8yj{k,T) 

Tr[Q2(fc)] 


Q,{k), 


(28) 


where Qj{k) were defined in (B 6 ). 

The functions yj{k,T) are solutions of the 
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(integro-)differential equations (for U = 0!) 


A. Scaling of the entanglement entropies 


drVjik^T) 


E 

e,nG{2,7,8} 



dp ye{p,T) 
2n el 


yn{k,T) 


+ h Y. (29) 

ne{2,7,8} 

where is the (3d) totally antisymmetric Levi-Civita 
symbol with = 1. The only nonzero elements of 
tensor A are given by 


Af{k) = 2 Af (A) 
Af(fc)=272 Afik) 


8 

1 + 7^ tan^ I 
87 ^ 

cot^ 1+72 


(30) 


while Qj read as 

5 ' 2 (fc) = 0, grik) = ■ gs{k) = l. (31) 

1 + ^ tan 2 

The initial conditions are as follows (c/. (26) and (28)) 


y2{k,0) = y7{k,0) = 0 

1 -cosfc 

ysik^O) = ±7 - - -. 


(32) 


where + and — correspond to /I = 0 and /3 = tt, respec¬ 
tively, in (25). 

The magnetic field appears as a multiplicative factor of 
the linear term in (29). This is because h is the coupling 
constant of , which is a quadratic operator in the J-W 
fermions (A12). Because only ys{k,0) is different from 
zero, the pre-relaxation limit is nontrivial only if /i 0 
(at the initial time the quadratic part of (29) is zero). 

We identify two non-interacting limits: A = 0 and 
A —> 0 at fixed Ah. The former is trivial (without per¬ 
turbation the expectation values can not move from the 
first plateau), while the latter corresponds to neglecting 
the quadratic part in (29). In the second case, (29) can be 
exactly solved and the solution is such that the Fourier 
coefficients of yj{k,t) relax to stationary values. This 
is the simplest nontrivial case of local relaxation in the 
pre-relaxation limit and was considered in Ref. [22] (see 
also the supplemental material [70]). Any displacement 
from the two noninteracting limits above is an effect of 
the interaction. 

It is important to note that this formalism can not be 
applied to the shift symmetrised state be¬ 

cause the mean-field description strongly relies on cluster 
decomposition properties. This is indeed a very clear ex¬ 
ample in which ignoring the importance of cluster decom¬ 
position has dramatic consequences: the (wrong!) initial 
conditions for the shift symmetrised state (i.e. the arith¬ 
metic mean of the initial conditions for (3 = 0 and (3 = tt) 
are identically zero, and hence time evolution would seem 
to be trivial in the pre-relaxation limit. Instead, consid¬ 
ering the equations for (3 = 0 ,tt separately, we will show 
that the shift symmetrizations of local operators like s| 
display nontrivial dynamics (cf. Fig. 2). 


We now focus on cases of local relaxation in the limit 
t ^ 0(A~^). Using the results of Ref. [69] we obtain the 
symbol of the mean-field Hamiltonian 


(k) = 


(I + m,(T))Q2(k) - jm,^s(T)Qsik) 

1+COSfe I 2 I —cosfe 

4 “*” / 4 


(33) 


mz(T) = l!v[p(T)(al -I- cr|)]/2 and m^^s(T) = 

Tr[p(T)(—tTf-|-cr|)]/2 are the magnetization and the stag¬ 
gered magnetization respectively 


m,{T)= f 

m^,s(T) = -7 / 

J — 


dfc y2(k,T) 


dk ys(k,T) 


27r 


St 


(34) 


One can verify that (k,T) generates the time evo¬ 
lution (29): idTT^^ik,T) = [n^^(k,T),r^^{k,T)]. By 
computing the characteristic length of the Fourier trans¬ 
form of (33) we infer®'^ that the typical range of is 
r^~ 2 /logl^l. 

Being essentially a noninteracting problem, one could 
wonder whether analytic expressions can be obtained for 
the time evolution of the entanglement entropies. There 
are two main complications: first, the solution of (29) is 
not analytically known and, second, computing the non¬ 
equilibrium time evolution of the entanglement entropies 
is not a simple task even in models that can be mapped 
to noninteracting fermions. The first goal is to decouple 
the two problems, which have different nature. 


1. On the explicit time dependence of 

The major obstacle to attacking the problem analyti¬ 
cally is that the mean-field Hamiltonian depends explic¬ 
itly on the time. In addition, the time dependence is not 
known a priori, being related to the solution of an infi¬ 
nite system of nonlinear differential equations. At first 
glance, this problem might appear insurmountable; in 
fact, it can be overcome as follows. 

We define the unitary operator 

Vrif) = e*^“^^C/MF(ir + T)C/^p(T), (35) 

where = H^^(oo) and 

z5rUMF(T) = H^^(T)Umf(T) . (36) 

In other words, Vt(T) generates mean-field time evolu¬ 
tion from T to T+T followed by a reverse evolution under 
^MF T+T to T . Let us assume that the mean-field 

Hamiltonian relaxes sufficiently fast for the existence of 
the following limits 

3 Jim Vt(T) = Vt 3 lim Vt(T) = I. (' 37 ') 

T->.oo T-foo > 
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We can then recast the problem into the non-equilibrium 
evolution under of the extrapolated density matrix 

m = VoP^^{0)V^. (38) 


instead justified only for the calculation of the extensive 
part of the entropies. 

After simple but tedious algebra we obtain the (ex¬ 
trapolated) correlation matrix r(fc, 0 ) of p(0) 


Indeed, using the existence of the limits in (37), we find 

= Vt[/mf(T) , (39) 

and hence 

p(T) = = Vtp^^(T)vJ: . (40) 

Let us analyze how Vr acts on operators. It is convenient 
to reinterpret Vt as a time evolution operator as follows: 

idfVrif) = S^T(f)VT{f). (41) 


f(/c,0) = hm 


tCtW 


8Xj(k) 

THQfM 

iG{2,7.8} 


E 


Q,(k) 


(44) 


with 


—2 sin lywz^s + i cos | (h -|- 2mz) 
^4sin^ + cos^ |(^ + 


(45) 


Here S^t(t) is a quasi-local effective Hamiltonian (with 
range increasing with r) that approaches zero in the limit 
of infinite r: 

jO^(r) = . (42) 

We denote by v(t) the (Lieb-Robinson) maximal veloc¬ 
ity associated with S)t(t) {v{t) is independent of T). At 
fixed time r, in an arbitrarily small time interval dr, the 
spreading of local operators under a quasi-local Hamil¬ 
tonian is proportional to the time interval and to the 
maximal velocity at which information propagates (cf. 
Ref. [91]). This means that Vt(T) is expected to map 
local operators into quasi-local operators with a typical 

diffusion length ^(T, T) ~ dru(r). We are going to 
assume that the limit ^(T) = limy_,^^ ^(T,T) exists and 
is hnite. This semiclassical estimate of ^(T) is a non¬ 
increasing function of T, so we find the upper bound 

aT) < m- 

Let us now consider subsystems with length £ ^(0). 

From (40) and our estimate of the diffusion length ^(0) 
it follows that the RDMs of p(T) are unitarily equivalent 
to the corresponding RDMs of p^^{T), up to boundary 
terms that extend over the typical length ^(0). This is 
independent both of the subsystem length and of T. On 
the other hand, the extensive part of the entanglement 
entropy is not affected by transformations that are uni¬ 
tary in the bulk, therefore, in the limit of large subsys¬ 
tem length, the entanglement entropies per unit length of 
must be equal to the ones of p{T). In conclusion, 
we mapped the problem into a sudden quench from p(0). 

We can sketch the steps of our reasoning as follows: 

(H- iTo)) ^ (H^F(T); pggEo) ^ p(0)), (43) 


and 


^ 2 )^) = Jim cos(2ffct) cosLlky2(k,i) 

i—>-oo 

- sinOfcCot ^ysikj)] + sin(2£fcf) , yrikj) 
2 J 2 sm ^ 


X 7 {k) = lim sm{2£kt) 

i—>-oo 


sin I 

—2cosOfc- y2{k,t) 


£k 


cos I 


2sinDfeii^ii^j/8(fc,0 + cos( 2 £kt)yr(k,t) 

£k 

k ~ 

Xs(k) = lim tan — sin Llky 2 (k, t) + cos ftkysik, t). 

i—^oc) 2 

(46) 

We called niz and mz^s the infinite time limits of rUzit) 
and mz^s(t) respectively; £k is the dispersion relation of 
the two species of quasiparticles diagonalizing 


£k = 


1 4 sin^ + cos^ I (^ + 2mz)2 

cos^ I -I- 7 ^ sin^ I 


(47) 


We note that the dispersion relation (47) becomes flat, 
and hence no more compatible with relaxation, when 
the magnetic field opposes the average magnetization re¬ 
stricted to even or odd sites h = 4z2mz^s — 2mz- 


2. Analytic expression for the entropy 

In noninteracting models the von Neumann entropy 
St = —Txpi log pe of i adjacent sites can be expressed in 
terms of the correlation matrix F^, restricted to the cor¬ 
responding subspace, as follows: 


where we used the notation (A; <I>) to indicate the dynam¬ 
ics generated by Hamiltonian A starting from state $. 
The first step is the mapping to a mean-field problem, 
which is justified by locality in the pre-relaxation limit. 
The second step is the reduction of the complicated 
mean-held dynamics to a sudden quench, which can be 


St = iTrJf'(f,), (48) 

with 

N 1 + a; , 1-1-a; 1 — a:, 1 — a: 

^(.x) = - — log ^- — log . (49) 
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The asymptotic behaviors of determinants and traces of 
(block-)Toeplitz matrices like have been extensively 
studied. We report here a strong limit theorem for the 
trace of a function / of an x iV block Toeplitz matrix 
T with smooth n x n symbol t{k) in the hypothesis that 
/ is analytic in an interval covering the spectrum of 

^Tr/(T) £ glTr/(t(fc)). (50) 


the interested reader to Section 3.1 of Ref. [9]) by the 
appropriate combinations of the 22,7,8; whose square is 
indeed proportional to the identity, like for the Pauli ma¬ 
trices. The only subtlety to take care of is that now the 
local space has dimension 4 and the length appearing in 
the expression is half of the subsystem length, being each 
point associated with two sites. 

For Tij ^ ~ At ~ 0(A°) we then find that the von 

Neumann entropy asymptotically behaves as follows 


This asymptotic expression can be used to compute the 
extensive part of the entanglement entropies when the 
symbol of the correlation matrix has no parameter com¬ 
parable with i. For 4-by-4 blocks we obtain 

< 51 ) 

We notice that this expression is invariant under uni¬ 
tary transformations on the symbol. Since the limit in 
(44) is assumed to exist, r(A:,0) is unitarily equivalent 
to F'^^(A:,0), so the effective sudden quench description 
gives in fact the same entropy (per unit length) of the 
GGE of the unperturbed model 

~Tr[Jf{T^^{k,0))] ■ ( 52 ) 

On the other hand, the scaling limit in which the 
rescaled time T = At is comparable with i can not be 
worked out using (50): the symbol depends on a param¬ 
eter comparable with the matrix size. This is the point 
where the mapping described in the previous section re¬ 
veals its importance: having mapped the problem into a 
sudden quench, we are in a position to apply the method 
proposed in Ref. [43] and slightly generalized in Ref. [9]. 
In particular, Refs [43] [9] considered (2£) x {2£) block 
Toeplitz matrices T with sufficiently regular 2-by-2 sym¬ 
bols t^^\k,t) of the form 

t^^\k,t) - +ni ■ , (53) 


s^t) 


AU 

-jr(Eo(*)) 


dk 


+ / — mm 

. 27r 


iin(4|4 




jr(Soo(fc))-^(So(fc)) 


(55) 


The Renyi entropies have the same func¬ 
tional form of (55), provided that the function is 

replaced with 

«^)-r^-((^)"-(^)")^ <-) 

The factor of 4 in the second line of (55) is in place 
of the factor of 2 in (54) as a result of a point being 
associated with two sites. One immediate consequence 
is that the actual maximal velocity at which information 
propagates (in the pre-relaxation limit) is double of the 
velocity associated with (47), i.e. vm = 2maxfc |5^|. The 
time at which linearity is lost is therefore given by 


tM = 7-7-• (57) 

4Amaxfe |t(.| 

The functions So(fc) and Soo(fc) are the eigenvalues of 
the symbol of the (extrapolated) correlation matrix at 
the initial and infinite (rescaled) times, respectively (no¬ 
tice that the 4-by-4 symbol of the correlation matrix has 
two doubly degenerate eigenvalues that differ only for the 
sign). In particular we have 


and worked out the limit I <C f ~ t for the trace of 
powers of T, obtaining the following result 


Tr[T2"(t)] 

2£ 


/ '^ dk f 

^-max(l-2|4|-,0)|nfcr 

/ '" dk f 

^-mm(I,2|e;|-)(n^)2", (54) 


where = {nlY + P- 

In our context translation invariance is realized by a 
two-site shift and, in turn, the correlation matrix has (at 
least) a 4-by-4 symbol. In fact this is not an issue be¬ 
cause, even in this case, the symbol follows an SU{2) 
dynamics (because of reflection symmetry), which is the 
main property used in the derivation of (54). The trick 
is essentially to replace the Pauli matrices that appear 
in the 2-by-2 symbol (53) of the original proof (we refer 


Eo(A) = 2, 


Ai(fc) , A?(fc) , Xlik) 


sin^ k 


sin" 


= 7 


sin I 

£k 


(58) 

where Xj{k) are given in (46) and the last equality reflects 
the equivalence with the unperturbed GGE (52). 

Despite Eq. (55) depending on dynamical parameters 
that are not known a priori, it is in fact an impressive 
result: the unknown variables do not affect the functional 
form of the entropy, which follows a linear behavior un¬ 
til tM and then saturates to a stationary value. This is 
exactly the type of universality that we were seeking in 
the late time evolution of the entanglement entropies. 

The physical interpretation depicted in Fig. 1 is per¬ 
fectly consistent with (55): Eg can be interpreted as the 
“cross section” for the production of a pair of compos¬ 
ite quasiparticles, consisting of two particles with almost 
the same velocity (the red and black arrows of Fig. 1, 
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treated as a single entity); Soo is instead related to the 
cross section associated with the constituents, which in 
the pre-relaxation limit are in fact distinguishable. 

Importantly, the necessity of replacing /9“^(0) by p{0) 
(38) can be interpreted as the effect of the scattering ma¬ 
trix not being trivial in the presence of interactions: we 
still describe time evolution in terms of free quasiparti¬ 
cles, but the interaction causes a dynamical modification 
of the effective cross sections entering into the problem. 

Finally, we propose a minimal refinement of (55) that 
combines the two scaling limits 1 ^ ~ t (computed in 

Ref. [43]) and 1 i At: 


Siit) 


^min(^4|4|^,l)^(Eo(/c)) 




dk 


in(4|4|^,l) [.^(So,(A:))-^(Eo(fc)) 


(59) 


This is the quantitative analogue of (6). 


3. A case study 


In order to use (55), we need another ingredient, 
namely the symbol of the (extrapolated) correlation ma¬ 
trix that emerges at infinite (rescaled) time. We will 
consider an example in which one-site shift invariance 
is restored, i.e. the staggered magnetization approaches 
zero ruz^s = 0 and the magnetization relaxes to a station¬ 
ary value. By inspecting (44) (45) (46) we find cosflk = 0, 
sinflfc = sgn(/i + ‘^niz), y 2 (k,t) becoming stationary and 
yr^sik,t) keeping oscillating with frequency 2Sk- 

First of all we evaluate the large time asymptotics of 
the maximal velocity v{t) of ^o(t) (42), so as to check 
whether the sudden quench description from the extrap¬ 
olated density matrix (38) can be used. The dispersion 
relation C£(fc,T) of i5o('r) is readily computed and reads 
as 


^ 2 / sin^ |7 ^w4(t) -f cos^ |(to4t) - 

y cos^ 1+7^ sin^ I 


(60) 

If Wz_,j(r) and uIz^t) — niz approach zero with different 
power laws, at long times we obtain 


v{t) ^ 2max(|TOz,s(r)|, |TO^(r) - m^j). (61) 


Thus, the typical diffusion length of Vr, namely 
40 ) ^ /g°°dru(T), is finite if the mean-field parameters 
relax faster than l/r. This will be the self-consistency 
condition to be checked a posteriori. 

Using that one-site shift invariance is restored, the 
symbol of the long time limit of the (extrapolated) corre¬ 
lation matrix is given by the one-site shift invariant part 
off(fc,0) (44): 


f{k, 00 ) = sgn(fe-|- 2 m^) Q 2 {k) 


sin k ej. 


2y2{k,oo) 
cos2 I el 


Q2{k). 

(62) 


From this we identify 

Eoo(fc) 


2y2{k,oo) 
cos |£fc 


(63) 


Remark Since the mean-field Hamiltonian depends ex¬ 
plicitly on the time, there are no simple shortcuts to 
compute the long time limit of dynamical quantities like 
(63). In particular Eoo(fc) can not be obtained from 
the correlation matrix of the GGE of the unperturbed 
model only knowing the mean-field Hamiltonian at late 
times. If we had approximated the dynamics by a sudden 
quench starting from the GGE of the unperturbed Hamil¬ 
tonian {i.e. r^^(fc, 0) instead of r(fc, 0)), we would have 
found the late time correlation matrix vanishing, namely 
Uoo(A:) = 0 and hence saturation of the entropies to the 
maximal values allowed®^. Such totally incoherent state 
oc I is equivalent to the GGE constructed with the one- 
site shift invariant conservation laws of the unperturbed 
model which in this case is also clearly equivalent to the 
GE. This is another example of the puzzling situation, 
pointed out in Ref. [69], where one-site shift invariance 
is restored in the pre-relaxation limit but the emerging 
stationary state is not constrained only by a deformation 
of the one-site shift invariant conservation laws of the 
unperturbed model. This could suggest that some extra 
degeneracy in the spectrum survives the leading order 
of perturbation theory. However, apparently this degen¬ 
eracy can not be associated with the non-abelian set of 
local conservation laws of the XY model, forming its one- 
site shift invariant charges an abelian subset. Ref. [69] 
proposed that this could signal issues in the limit A —>■ 0 
of quasi-local (quasi-)conserved operators, but, in fact, 
this is still an open question. 


B. Numerical analysis within mean-field 

1. Prediction vs numerical solution. 

Formula (55) is written in terms of quantities com¬ 
puted at late times after the quench, which must be ex¬ 
tracted from the numerical solution of the equations. 

In Fig 2 we show the time evolution of the magneti¬ 
zation, its time derivative and the staggered magnetiza¬ 
tion in the pre-relaxation limit for a particular quench 
with 7 = 0.75, h = 0.75 and U = 0. The numerical 
data leave no doubt the parameters of the mean-field 
Hamiltonian time relax and one-site shift invariance is 
finally restored. In particular, our best estimate of the 
late time magnetization is rriz ~ —0.0672556. We notice 
that both mean-field parameters relax faster than 1/At, 
corroborating the validity of the mapping (43) to a sud¬ 
den quench from the extrapolated density matrix (38). 
The asymptotic value of rriz determines the time tM (57) 
at which the entire subsystem becomes (substantially) 
entangled with the rest, but the full time evolution of 
the entanglement entropies requires also the knowledge 
of Eoo(fc). 
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FIG. 2. (Left)For the quench (a) displayed in the panel, the magnetization in the 2 direction relaxes to a stationary value. The 
inset shows the time derivative of the magnetization approaching zero as (Right) The staggered magnetization decays 

as (At)“5. For the sake of clarity, in the plots in log-log scale the regions under the curves are shadowed. 



FIG. 3. The quantity EAt(fc) = as a function of 

cos 

the momentum k for 4 times after the quench (a) displayed 
in the panel. The curves become thinner as the time gets 
larger. Gonsistently with (46) when one-site shift invariance 
is restored, the limit of infinite time for y 2 {k,At) seems to 
exist. 


Fig. (3) shows the quantity T^rik) = as a func- 

COS 2^k 

tion of the momentum k for various times*^^. The limit 
limT-j-oo ^^(A:) = Soo(A:) is reached rather quickly, at 
least for momenta not too close to ±7r. 

Having checked that the parameters 7 = 0.75 and 
h = 0.75 are compatible with relaxation and restoration 
of one-site shift invariance, we can now compare the mean 
field solution with the asymptotic formula (55). In Fig. 4 
the time evolution of the entanglement entropy per unit 
length as a function of the rescaled time At is plotted for 


various subsystem lengths. At any given rescaled time, 
the extrapolation is done by finding the best parabolic 
fit in l/£, weighted by (the inverse of the square of the 
leading correction neglected). 

The excellent agreement between the left panel of 
Fig. 4 and (55) confirms the conjecture that Vr, (40), 
is only responsible for a quasi-local mixing that does not 
affect the extensive part of the entanglement entropy. We 
also notice that the value of the entropy per unit length 
associated with the second plateau is visibly smaller than 
log 2 , which would have been obtained neglecting the 
mean-field dynamics at intermediate rescaled times (be. 
approximating Vr, (37), by the identity). 

Figure 4 exhibits the scaling behavior claimed in ( 6 ). 


2. Intermediate mean-field 

The results of the previous sections, culminating in 
(55) and Fig. 4, rely on the mean-field description, which 
is supposed to emerge in the limit A <C 1 at times 
t ^ A“^. In the pre-relaxation limit it was legitimate to 
replace the initial state by the GGE of the unperturbed 
model, however in actual numerical simulations this in¬ 
troduces relevant corrections that can spoil the agree¬ 
ment with the asymptotic values. Here we discuss the 
extent of such corrections by considering an intermedi¬ 
ate mean-field approximation, where the initial state is 
not replaced by p^^(O), as instead done in (A9). We 
consider the following system: 

(C’)mf (0 = (^mf(A)|C>|4'mf(A)) 

tdt |4/MF(t)) = [Ho + AH^^iAt)] |4'mf(A)) (64) 

|4'mf( 0)) = Idio) • 
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FIG. 4. The von Neumann entropy per unit length as a function of the rescaled time per unit length for various intervals 
within mean-held. At any given rescaled time, the extrapolation is done by hnding the best parabolic ht of the data, weighted 
by (being too close to the asymptotic values, the data for ^ = 128 are not displayed in the left panel). The corresponding 
conhdence interval at 95% is shown as an orange region (in the left panel the width of the region has been increased to make it 
visible). (Left) For the quench (a) the agreement with formula (55) is excellent. The gray dashed vertical line is (57). (Right) 
Another quench (b) that leads to local relaxation (in the pre-relaxation limit) though one-site shift invariance is not restored 
and there are more pronounced oscillations; data are still compatible with the scaling law (6). 




FIG. 5. The time evolution of the von Neumann entropy per unit length within the intermediate mean-field (Gl) as a function 
of the rescaled time per unit length for the quench (a) displayed in the panels and various values of A. The black curves 
corresponding to A —>■ 0 are displayed in the left panel of Fig. 4 in purple (left) and red (right). The gray vertical line is (57). 
In the left panel the dashed lines are the refined predictions (59). 


Differently from (A9), we have not written the equations 
in terms of the rescaled time T = At. This choice re¬ 
flects the fact that the mean-field time evolution of j'l'o) 
includes contributions that vary in a timescale ^ A*^. 
The immediate advantage is that such description is prac¬ 
tically undistinguishable from exact time evolution at 
times t <?; A“^, when the Hamiltonian can be approxi¬ 
mated by Hq. Then, in the pre-relaxation limit, it be¬ 


comes equivalent to (A9)®®. We do not expect any im¬ 
provement in the corrections o(A°) at times t ~ A“^, 
which are in fact beyond mean-field, but the interme¬ 
diate mean-field is useful to recognize the effects of the 
preliminary relaxation to the first plateau. Incidentally, 
we point out that we worked out system (A9) instead 
of (Cl) because in the latter the number of differential 
equations is doubled and half of the information is useless 
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FIG. 6. Evolution of the average magnetization rriz after the quenches (a) and (b) displayed in the panels for various “small” 
values of A. The symbols are the iTEBD numerical results and the thin colored lines correspond to mean-field (ME) time 
evolution (Cl) starting from the initial state. The thick red lines are instead the ME predictions (A9) starting form the 
unperturbed GGE (dotted lines). For the sake of clarity, in the left (right) panel both datasets and MF curves corresponding 
to A > 0.05 (A > 0.025) are vertically shifted (the red curves are in fact coincident!). 


in the pre-relaxation limit, making the solution of (Cl) 
less indicated for analytic investigations. In addition, we 
are able to compute the time evolution (Cl) only if the 
initial state is a Slater determinant, whereas in princi¬ 
ple (A9) can be used also quenching from an interacting 
model. 

The left panel of Fig. 5 shows the time evolution of the 
entanglement entropy per unit length under (Cl) for var¬ 
ious values of A and fixed subsystem length £ = 32. The 
scaling behavior of (6), including the preliminary linear 
growth in the space-time scaling limit 1 i < t/{2vM), 
is completely revealed. In addition, the data are in good 
agreement with the prediction (59). It is however impor¬ 
tant to note that the actual corrections to the mean-field 
solution could spoil the physical meaning of the plot. 
The larger A is and the shorter the time window that 
can be described by mean-field: it is unreasonable to ex¬ 
pect agreement with exact dynamics when At ~ A“^! 
For the subsystem length considered in the right panel 
of Fig. 5 (t = 4) these potential problems should be less 
severe. Nevertheless, oscillations apart, the behavior is 
essentially the same as in the left panel. 

C. Mean field vs iTEBD simulations 

This section is devoted to comparing the mean-field 
predictions with the exact quench dynamics induced by 
Hamiltonian (Al) and obtained via iTEBD numerical 
simulations for quenches starting from the Neel state 
(some details of the numerical simulations are reported 
in the supplemental material [70]). We consider two dif¬ 
ferent quenches, with parameters: 


(a) 7 = 0.75, h = 0.75, U = 0; 

(b) 7 = 2, h = 1, H = 0.25. 

In both cases we expect relaxation in the pre-relaxation 
limit but in (b) one-site shift invariance is not restored. 
In the iTEBD simulations the (small) parameter A has 
been chosen in the interval [0.0125,0.2]. The computa¬ 
tional complexity does not change much with the interac¬ 
tion strength, making it extremely difficult to reach large 
rescaled times At for the smallest values of A. This is 
the reason why we have also examined quenches with val¬ 
ues of the interaction that can not be considered strictly 
small, e.g. A = 0.1 or 0.2. Nonetheless, we would empha¬ 
size from the very beginning that, even for these values of 
A, we expect an overall trend that should at least qual¬ 
itatively agree with mean field. In other words, even if 
the crossover is mathematically defined only in the pre¬ 
relaxation limit, the effects should remain visible also for 
finite values of A. 

Fig 6 shows the time evolution of the magnetization 
rriz. The iTEBD numerical data (symbols in the plots) 
exhibit deviations from mean-field that are reduced as A 
gets closer and closer to zero. This is exactly what we 
expect: a finite value of the interaction strength neces¬ 
sarily introduces some corrections to mean-field dynam¬ 
ics. The magnetization moves away from the GGEq value 
rriz = 0 and, up to subleading corrections, data for dif¬ 
ferent A and fixed At collapse to a curve that is compat¬ 
ible with the mean-field prediction. This is a conclusive 
confirmation that the pre-relaxation limit is nontrivial 
and apparently described by mean-field. In addition, the 
magnetization seems to approach a nonzero value for ar¬ 
bitrarily small A. Since the energy density is — ^ and, in 
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FIG. 7. Evolution of the von Neumann entropy for a subsystem with 1=1 (top) and 1 = 2 (bottom) sites after the quenches 
(a) and (b) displayed in the panels for various “small” values of A. The notations are the same as in Fig. 6. 


turn, in the limit A —>■ 0 the thermal expectation value 
approaches zero, one can infer just from numerics that 
the second plateau emerging in the pre-relaxation limit 
is different from the infinite time limit. It corresponds to 
the highly nontrivial pre-thermalization plateau that, as 
mentioned in the introduction, is generally developed af¬ 
ter a pre-relaxation plateau when the unperturbed model 
is non-abelian integrable. 

In Figs 7 and 8 we show the time evolution of the en¬ 
tanglement entropy for subsystems with £ = 1,2, 4, and 
6 lattice sites, respectively. In particular we compare the 
iTEBD data with both the mean-field prediction starting 
from the GGEq (thick red lines in the figures) and the 
intermediate mean-field (Gl) (thin colored lines). Let us 
analyze more attentively the behavior of the entropies. 
Generally, S(, experiences an almost linear growth (apart 
from oscillations) up to a time t* ~ £/{2vm), where the 
maximal velocity is given by um = |7 ~ 1| + o(A°); af¬ 
ter that, the entropy approximately reaches the corre¬ 
sponding GGEq value. This transient represents the well- 


known Galabrese-Cardy regime^ of linearity after a global 
quench. Obviously, in the limit A —)■ 0 at fixed rescaled 
time At, this is squeezed into an infinitesimal slice at zero 
rescaled time (c/. Eig. 5), so the entropy in the GGEq 
becomes the starting value of the entropy in the mean- 
field description. However, the iTEBD simulations are 
done at fixed A and the larger the interaction and the 
subsystem size are and the more evident the transient is. 

We observe a good agreement between numerical sim¬ 
ulations and mean-field, within the limits of corrections 
that tend to zero as A —)■ 0. Nevertheless, some com¬ 
ments are in order: the interplay between i and A de¬ 
termines the size of the corrections to mean-field. We 
expect mean-field to be valid in a regime At ^ A““ and 
£ A~^, with a and /3 some unknown exponents. The 
results of Ref. [68] suggest that thermalization in generic 
models is expected in a timescale ~ A“^, so we can set 
a = 1; concerning /3, the emergent mean-held scaling 
regime At ~ £ points to a correspondence between £ and 
At, so we assume /3 ~ a = 1. These estimations allow 
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FIG. 8. The same as in Figure 7 for a subsystem with £ = 4 (top) and £ = 6 (bottom) sites. 


us to understand why for Sq the data show larger de¬ 
viations from mean-field: for A = 0.2 (or 0.1), one has 
A£ ~ 1, which could be outside the region of validity of 
mean-field. On the other hand, Si and S 2 show excellent 
agreement with the theoretical predictions. 

In order to establish the validity of the mean-field dy¬ 
namics we also carried out a meticulous analysis of the 
scaling properties of the normalized Frobenius distance 

'DiPe^pT^) = \\pe - pf^\\/{\\pe\? + Wpf^W^Y^^ consid¬ 
ered in Ref. [11] (see also the supplemental material [70]) 
between the iTEBD reduced density matrices and those 
obtained via (Cl). 

In Fig 9 we show the time evolution of T>{pi,p^^) for 
different values of i and A. In the regime of validity of 
mean-field A£/(2|7 — 1|) < At < A~^ (the lower bound 
is for the mean-field (A9)), the distance I?(pf,p^^) is re¬ 
duced in magnitude as A becomes smaller and smaller. 
Moreover, for fixed A, the distance increases as the sub¬ 
system length becomes larger. 


If the assumptions behind (A9) are satisfied, for A —>• 0 
the distance between the two reduced density matrices 
should vanish. Showing the base-A logarithm of the dis¬ 
tance for a sequence of A at fixed ratio between the one 
and the next (A = 0.2 • 2“") is the way that we chose to 
bring it out: a nonzero limiting distance would result in 
the logarithm to approach zero as ~ l/(n -I- 2.3), which 
can be fairly recognized by sight (instead, the logarithm 
remains nonzero if the distance scales as a power law in 
A). Our impression based on the analysis of the insets 
of Fig. 9 is that the data are more compatible with the 
distance approaching zero as a power law in A instead of 
reaching a nonzero value that might have signaled some 
issues in the hypotheses behind (A9). 


VII. PERSISTENT OSCILLATIONS 

As briefly discussed in Section IV B, there are situ¬ 
ations in which the solution of the mean-field equations 
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FIG. 9. (Main) The normalized Frobenius distance (see the supplemental material [70]) between the exact RDM and the one 
obtained via the intermediate MF description (Cl) for subsystems with a number of sites ^ = 2, 4, 6 for the quenches (a) and 
(b) displayed in the panels and various “small” values of A. As A gets smaller, at fixed At the distance is reduced. (Inset) 
The logarithm of the distance in base A is shown. The data suggest the logarithm is approaching a finite value, that is to say 
the distance decays to zero as a power-law in A. 


displays persistent oscillatory behavior at arbitrarily long 
times. In Ref. [69] persistent oscillations have been inter¬ 
preted as manifestations of localized (quasi-)excitations. 
Following that point of view, the oscillatory behavior is 
not expected to modify the extensive part of the entangle¬ 
ment entropy. We test this conjecture on two quenches: 

(c) Time evolution of the ground state of the Majumdar- 
Ghosh Hamiltonian (/3 = ^ and ^ = tt in (A4)) under 
the Hamiltonian (Al) with 7 = 0.5, h = 0.5 and 
U = 0.5. 

(d) Initial state that breaks reflection symmetry, namely 
(3 = ^ and (/> = tt in (A4), with Hamiltonian param¬ 
eter 7 = 2, h = 0 and U = —1.5. 

We briefly comment on the mean-field solution of the dy¬ 
namics. In case (c) the nontrivial part of the dynamics 
involves the families of charges with symbols (B 6 ) propor¬ 
tional to Q 2 {k), Qs]^), and Qeik). Instead, in case (d) 
the relevant symbols are Q 4 {k), Q^{k), and Qs{k). To¬ 
gether with the quenches from the Neel state (where the 
symbols involved were Q 2 {k), Qrik) and Qs{kj), these 
are representative of all the quenches starting from states 


of the form (A4) in which the mean-field dynamics is re¬ 
duced to three families of local conservation laws. We 
do not report the corresponding systems of differential 
equations, which however can be found in Ref. [69]. 


A. Numerical analysis within mean-field 

Figure 10 shows the time evolution of the von Neu¬ 
mann entropy for the quenches (c) and (d) . In both cases 
the oscillations enter as subleading corrections to the en¬ 
tropies per unit length, as it was conjectured in Section 
IV B. In the left panel the nontriviality of the scaling limit 
is manifest. In addition, the data are compatible with the 
asymptotic formula ( 6 ). We notice however that the in¬ 
crease in the entropy is modest. This is not unexpected 
since the presence of persistent oscillations signals that 
the dephasing mechanisms behind the relaxation process 
are less effective. As a consequence, the entropy growth 
is restrained. This effect can be so important that the 
extensive part of the entropy could appear almost inde¬ 
pendent of time, as in the right panel of Figure 10. In 
the limit of large subsystem length the entropy practi- 
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FIG. 10. Evolution of the entanglement entropy per unit length within mean-field (A9) for various lengths for the quenches 
(c) and (d) displayed in the panels. In both cases, expectation values of local observables show persistent oscillatory behavior. 
(Left) In the pre-relaxation limit the entropy per unit length exhibits the same scaling behavior (6) as in the cases of local 
relaxation. (Right) Persistent oscillations are an indication of less effective dephasing mechanisms: for large subsystems the 
entropy per unit length does not grow significantly. The red curve (low-pass) is an extrapolation based on the average of the 
upper and lower envelopes. The accuracy of the extrapolations is clearly spoiled by the sizable oscillations. 


cally does not move from the value associated with the 
unperturbed GGE. This time, even though very large 
subsystems are considered, it is difficult to identify the 
linear part of the time evolution. However, an extrapola¬ 
tion with confidence level 95% suggests that the scaling 
limit is likely to be different from the value corresponding 
to the unperturbed GGE. 

Since we have not analytically computed the asymp¬ 
totics of the entropy for quenches with persistent oscil¬ 
lations and our intuition is only based on the qualitative 
picture of Section IV, we can not rule out the emergence 
of other scaling behaviors when there are persistent oscil¬ 
lations and the initial state is sufficiently general. Specif¬ 
ically, the mapping of the mean-field time evolution into 
a sudden quench described in Section VIA1 can not be 
used without modifications and, for example, we can not 
exclude that the limit At —>■ 0 in the asymptotic expres¬ 
sion could be different from the value corresponding to 
the GGE of the unperturbed model. The right panel of 
Fig. 10 suggests indeed this might be happening. Never¬ 
theless, all the cases that we considered are compatible 
with the scaling law (6). 


B. Mean field vs iTEBD simulations 

In this section we show the results of the iTEBD sim¬ 
ulations for the quenches (c) and (d). In both cases we 
focus our attention on the entanglement entropies for 
£ = 1, 2, 4, and 6 and on the distances between the re¬ 
duced density matrices of the exact dynamics and the 
reduced density matrices obtained through the interme¬ 


diate mean-field approach. Below we list the results of 

such comparison together with some considerations: 

(c) The Majumdar-Ghosh Hamiltonian ground-state 
[vj/g) = |MG) = {^(Iti) ~ l'l't))/'\/2 breaks transla¬ 
tional invariance. Notwithstanding, the local magne¬ 
tization (cr|) /2 is translational invariant at any time 
after the quench. This is a consequence of reflection 
symmetry abont a bond. Therefore, we do not show 
the time evolution of the magnetization, as it gives no 
more information than the single-site entanglement 
entropy. 

In Fig 11 we show the evolution of the von Neumann 
entropies. The agreement with the mean-field pre¬ 
diction is pretty good and improves as the system 
length and the interaction strength get smaller. The 
discrepancy is more evident at larger rescaled times 
and could signal the emergence of another regime for 
At ~ A“^. Essentially, the same considerations ap¬ 
ply to the distance between reduced density matrices 
(see Fig 13). This remains rather small (especially for 
£ = 2) even for the largest value of the interactions we 
considered (A = 0.2). Finally, also in this case the 
hypothesis of a power-law scaling for the distances 
seems to be confirmed by the plots of the logarithm 
of the distance in base A. 

(d) The initial state ITq) = 0(-\/3|ti) ~ lit))/2 breaks 
both one-site shift invariance and reflection symme¬ 
try. We do not show the average magnetization 
which for h = 0 is zero at any time. This is a con¬ 
sequence of the combined symmetry under reflection 
about a bond followed by spin flip. 
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FIG. 11. Evolution of the von Neumann entropies for the quench (c) displayed in the panels. For the sake of clarity all datasets 
with A > 0.025 have been vertically shifted. The notations are the same as in Fig. 7. 


The evolution of the von Neumann entropy is de¬ 
picted in Fig 12. It is evident that the corrections 
to the mean-field approximation are stronger than 
in the previous cases; nevertheless, they still reduce 
as A gets closer to zero. Interestingly, due to the 
absence of a magnetic field, the approximation of ne¬ 
glecting the interacting terms for this specific quench 
is trivial, reducing to no evolution at all. There¬ 
fore, the nontrivial dynamics of the iTEBD entropies, 
which undoubtedly move away from their initial val¬ 
ues, show unequivocally that the pre-relaxation limit 
does exist and is not trivial even for a genuine four- 
fermion perturbation. 

The overall behavior and the scaling with respect to 
^ and A is confirmed by the analysis of the distances 
between the reduced density matrices (see Fig 13), 
which basically follow the same qualitative behavior 
of the difference between the iTEBD entropies and 
the mean-field ones. 


VIII. CONCLUSIONS 


In this work we have considered the evolution of the 
entanglement entropies after a quantum quench from ini¬ 
tial states that break one-site shift invariance in a per¬ 
turbed quantum XY model. Our main result is to have 
identified a long-time behavior that is almost indepen¬ 
dent of the quench details. We showed that, in the 
limit t ~ A~^, where A is the strength of the pertur¬ 
bation, the entropy, which already reached a plateau 
of prethermalization/pre-relaxation, grows linearly in At 
until a characteristic time Im ~ after which it starts 
saturating to another quasi-stationary value. We pro¬ 
posed a semi-classical explanation, based on very gen¬ 
eral arguments, which manifests the universality of the 
phenomenon. The nontrivial behavior disclosed in the 
timescale At ~ 0(A°) is driven by the breaking of what 
we called “non-abelian integrability”, namely the exis¬ 
tence of an infinite set of local conservation laws that do 
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FIG. 12. The same as in Figure 11 for the quench (d) displayed in the panels. 


not commute with one another. When the perturbation 
removes the accidental degeneracies responsible for the 
non-commuting local conservation laws, at the leading 
order in A for t ~ local observables become sensible 
to the slow unitary time evolution that occurs in each 
(quasi-)degenerate subspace. As one would expect from 
general reasonings, the entanglement entropy reacts to 
the symmetry breaking by increasing its value. 

We have studied in details a quench from the Neel 
state under the Hamiltonian of the XYZ spin-^ chain 
with an external magnetic field that breaks integrabil- 
ity. In particular, we obtained an analytic expression 
for the time evolution of the entanglement entropies per 
unit length in the scaling limit 1 <C £ ~ At. Considering 
other quenches, we showed that the pre-relaxation of the 
entropy is almost independent of the quench parameters 
and, even when local observables exhibit persistent os¬ 
cillations, in the entanglement entropies per unit length 
the oscillations turn out to be subleading. 

Our analytic findings are based on the method devel¬ 


oped in Ref. [69] : the dynamics are reduced to mean-field 
time evolution and the entropy is computed using free- 
fermion techniques. 

We checked our results against iTEBD simulations. 
We had to face the technical problems behind the pre¬ 
relaxation limit, which presupposes very long times. The 
present-day algorithms are reliable on timescales that, in 
the cases considered, do not depend much on the strength 
of the perturbation, so the smaller A is and the shorter 
the maximal rescaled time At that is reached. In order 
to make the comparison more transparent, we analyzed a 
distance between the reduced density matrices that time 
evolve with the exact Hamiltonian and those of the mean- 
field description. Notwithstanding the complications of 
the limit, we provided some numerical evidence that the 
mean-field time evolution is exact in the pre-relaxation 
limit, confirming that the technical assumptions made in 
Ref. [69] are reasonably satisfied for the time evolution 
under Hamiltonian (Al). 

Our work raises a number of issues. First and foremost 
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FIG. 13. The normalized Frobenius distance (see the supplemental material [70]) between the exact RDM and the one obtained 
via the intermediate mean-field (Cl) for the quenches (c) and (d) displayed in the panels. The notations are the same as in 
Fig. 9. 


is the perceptible inadequacy of DMRG-like algorithms 
to study time evolution in the pre-relaxation limit or in 
other limits in which the time is supposed to be very 
large. Alternative methods based on a perturbative so¬ 
lution of the equations of motion have been recently em¬ 
ployed to study similar problems®®. Apparently these 
techniques allow one to reach larger times, although the 
errors are much less under control than in DMRG algo¬ 
rithms. It could be very interesting to investigate the pre¬ 
relaxation limit with these methods, checking the consis¬ 
tency with mean-field. 

Another relevant question is whether the time- 
dependent GGE (18) that emerges in the pre-relaxation 
limit survives the next order of perturbation theory, 
like a (deformed) GGE description does at intermedi¬ 
ate times in the presence of small integrability breaking 
perturbations^'. 

Finally, our analysis suggests that different local con¬ 
servation laws can respond to perturbations in different 
timescales. The pre-relaxation time Im seems to char¬ 
acterize the time at which some non-commuting con¬ 
servation laws of the unperturbed model with range I 
start feeling the effects of the perturbation. Other local 
charges are instead expected to remain quasi-conserved 
for much longer times. A systematic study of these ef¬ 


fects could be useful to identify the correct set of con¬ 
servation laws that contribute to the description of pre- 
relaxation/prethermalization plateaux. 

We conclude with a speculation. Our findings rely on 
the non-abelian integrability of the unperturbed model, 
which in our case is noninteracting. A natural question 
is if there are also interacting quantum many-body sys¬ 
tems with infinitely many non-commuting (quasi-)local 
conservation laws. Although there is not yet firm ev¬ 
idence, models with a loop algebra symmetry, like the 
XXZ model with anisotropy root of unity®®, are good 
candidates for possessing similar charges. This is clearly 
an imperative question to be addressed in the next fu¬ 
ture. 
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SUPPLEMENTAL MATERIAL 
Appendix A: On the iTEBD simulations 

In this supplementary section we discuss the numerical method used to access the pre-relaxation limit of the exact 
quench dynamics under Hamiltonian 


H = Ha + A{US^-^ + hS ^), 


(Al) 


with 


Ha = 


E 

t 


^ ^ -X X I 1 y y , ^ z z 
^ ^£+1 ‘^£‘^£+1 + ^^^£^£+1 ) 


(A2) 




ji: 




1 

2 



(A3) 


and initial state 


l^-o) =' 


cos^ Iti) -he^'^sin^ lit) 


(A4) 


In particular, we make use of intensive numerical simulations via infinite time-evolving block-decimation (iTEBD) 
algorithm®®. 

The algorithm is based on the matrix product state (MPS) description of one-dimensional lattice models and works 
directly in the thermodynamic limit. Since the Hamiltonian (Al) includes interactions between next-nearest neighbor 
spins, we implemented the algorithm by considering a unit cell of two sites®^, with local dimension d = 4. The initial 
state has then a translational invariant MPS representation. Nevertheless we need to partially break translational 
symmetry in order to simulate the action of the local evolution operators on even and odd bonds. A generic many-body 
state takes the following MPS canonical representation 


Id/) = ^ Tr[... A„r^/+/Ae- • • ] I- • • • •). 

IH 


(A5) 


Here are x ^ X matrices associated with odd/even unit cells, with Sj spanning the j*® unit-cell space with 
canonical basis {|tt) i Iti) i lit) j lit)}; similarly, A^/g are diagonal matrices with entries the singular values associated 
with the rank-x Schmidt decomposition for the bipartition of the system on the odd/even bonds. The trace is taken 
over the auxiliary space. We point out that (A5) is well-defined thanks to the right/left orthonormalization conditions 


^(r:/gA,/g)(r:/gA,/g)t = I, ^(Ag/,r:/g)t(Ag/,r:/g) = i. 


(A6) 


which essentially state that the leading eigenvalue of the right/left transfer matrices is unitary and the corresponding 
right/left eigenvector is the identity matrix (reshuffled as a vector). These conditions guarantee the correct normal¬ 
ization of the vector Id/) as well as the possibility to perform well-defined operations in the MPS representation. 

The initial state |d/o) (A4) admits a trivial rank-1 MPS representation with a unique singular values = 1 and 
a 1-by-l tensor = cos f sin f 

The time evolution of the MPS in Eq. (A5) is easily accomplished by using the second order Suzuki-Trotter 
decomposition of the evolution operator with time step dt, namely 


^—iHdt 





j odd j even j odd 


(A7) 


where the index j runs over the unit cells and t) represents the local interaction between nearest neighbor cells. Since 
the original Hamiltonian is translational invariant, the local interaction 1) does not depend on the particular bond 
between cells. In practice, the MPS representation of the state at time t + dt is obtained by updating the MPS at 
time t so as to take account of the repeated action of the local evolution operators in Eq. (A7). Interestingly, under 
the action of the unitary operator restricted to the odd (even) bonds, only the matrices {Pg,Pg, Ag} ({Pg,Pg, Ag}) 
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need to be updated and the invariance of the evolved state under a shift by two unit cells (z.e. four chain sites) is 
preserved. 

In our numerical simulations we fix dt = 0.05 (we verified that the data are not affected by the time discretization). 
The auxiliary dimension y used to describe the reduced Hilbert space is dynamically updated in such a way that, at 
each local step, all the Schmidt vectors corresponding to singular values larger than Amin = 10“^® are retained. For 
purely practical reasons, this condition is relaxed when y reaches a maximal value Xmax G [400,800] dependent on the 
particular quench. Because of the (unavoidable) upper bound Xmax, the truncation procedure is the main source of 
errors of the algorithm. We stop the simulations so as to have a Schmidt error coming from the iterative truncation 
of the Hilbert space always smaller than 10“^. We point out that the bipartite entanglement entropy grows linearly 
for all explored times, as it should do after a global quantum quench. 

Based on the normalization condition (A6), the evaluation of expectation values of local quantities with spacial 
range £ turns out to be very efficient and scales as £x^- However, for our purpose, we need to reconstruct the full 
reduced density matrix pg for a subsystem with £ adjacent sites, which is actually computationally more demanding. 
In particular, for £ even one has 


Pi 


Tr[(Aer^ 


o/e 0/6 


)^(Aer: 


■r;/;A„/,) 


(A8) 


where the last two matrices appearing in the strings can be odd or even depending on the parity of £/2 and the trace 
is over the auxiliary space. Clearly, the reduced density matrix of an odd number of sites is obtained by tracing out 
the last spin, which corresponds to trace (A8) over the second site in the last unit cell. Finally, from the spectrum of 
the reduced density matrix one can easily evaluate the von Neumann entropy Si = —Tr(p£ logpi) or, more generally, 
the Renyi entropies 5'^“^ = log[Tr(p“)]/(l — a). 


1. Distance between reduced density matrices 

The main difficulty in comparing the theoretical predictions with the numerical data is that, in order to describe 
the scaling limit A —0 with At fixed, the simulations should reach very large times t ^ 1. In practice, all algorithms 
based on MPS representations suffer from the growth of the bipartite entanglement entropy, making it extremely 
challenging the possibility to investigate the pre-relaxation limit. This forces us to choose not too small values of 
the parameter A, introducing unavoidable corrections o(A*^) which are beyond our control. The agreement with the 
mean-field approximation is therefore made evident via a systematic scaling analysis. As A becomes smaller and 
smaller we aim at showing that, for fixed value of the rescaled time At, the expectation values of local quantities 
converge to the mean field predictions; more generally, the full reduced density matrix pi should get closer and closer 
to the corresponding reduced density matrix p^^ of the mean-field state defined as follows 

(^^)gL„ (T)=Tr[p“^(T)0] 

^^TP^^(T) = [H^^{T),p^^{T)] (A9) 

p^^ (0) = lim lim Tr^ [e*^“‘ | (4-01 . 

|>l|->.oo t->oo 

The approach to the mean-field solution strongly depends on the subsystem size £. Indeed, in the mean-field 
approximation the condition 1/A ^ is always verified (since A has been already sent to zero). On the other hand, 
in the iTEBD numerical simulations, for a finite small value of A, one needs to consider small enough subsystems to 
be in the expected regime of validity of mean-field. 

In order to give a quantitative estimation of the deviations between the numerical simulations and the mean-field 
results it is convenient to compute a distance between the corresponding reduced density matrices. Ref. [11] analyzed 
the physical meaning of different distances and found a good compromise between practicality and meaningfulness in 
the following normalized distance: 


{P: P) 


IIp-pII 
y/\W + \\p\? 



2Tr[pp] 
Tr[p2] + Tr[p2] 


(AlO) 


where jjOjj = \/Tt(0^0) is the Frobenius norm of an operator. The distance (AlO) has several useful properties'^; 

(i) It vanishes if and only if the reduced density matrices are identical. 

(ii) It takes values in the interval [0,1]. 
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(iii) It has an interesting lower bound T>{p,p') > \e ^ s *'^V 2 |^^g-S( 2 ) _j_g-S'( 2 )^ 

(iv) Even in the limit of large subsystem, if the density matrices are physically different, the distance does not 
approach zero (this is a problem that affects other distances, like the operator norm^^). 

Notice that Property (iii) implies that it is sufficient that the Renyi entropies of the two RDMs differ by a 
constant for V{p,p') to have a lower bound independent of 

Since the RDMs of the exact dynamics are not gaussian, one can not use free-fermion techniques to compute (AlO) 
but instead the RDMs within mean-field must be explicitly reconstructed. In particular we have 



mf 

Tiiexp{ajKff^ai)] ’ 


(All) 

where aj are the Majorana fermions 




02^-1 = n CTiO-? a2i = 

j<i j<e 


(A12) 

and the (implicit) sums are 
the correlation matrix TN^^ 

•'J 

restricted to i adjacent sites; the matrix elements are related to the elements of 

= 5ij — {aiaj)MF via the matrix relation = arctanh(r^^)/2, where T^^ is the 


correlation matrix restricted to £ sites. These properties are a consequence of the gaussianity of the many-body state 
in the mean-field description. 


Appendix B: Pre-relaxation in quantum XY models 


In this section we integrate the main results reported in the paper with the pre-relaxation limit of the entanglement 
entropies in the noninteracting case described by the XY Hamiltonian Hq with a small transverse field A/i 


^ = E 


1 -t- 7 


1-7 

4 


y y I A/i 


(Bl) 


The time evolution under (Bl) of a two-site shift invariant Slater determinant was worked out in Ref. [22]. We 
reconsider the same problem here. 

It is known^^ that in noninteracting models the symbol r(fc) of the correlation matrix time evolves with the symbol 
7i{k) of the Hamiltonian as follows 

r(fc, t) = . (B2) 

In our specific case H(fc) = 'Ho{k) + Ahl (g) so we have 

T{k, t) = V\k, t)e-*’^Afe)tr(/!c, V(A;, t) (B3) 

with 

idtV{k, t) = -A/ie-*’^“g V{k, t). (B4) 

The mean-field approximation Yupik, ht) ^ V{k, t) corresponds to replacing by its time average, 

that is to say (for 7 ^ 0 , otherwise the perturbation commutes with Hq) 


I4if(^, ht) = e 

with Q 2 {k) the symbol associated with the second family of local conservation laws of Hqj namely 

Qi(fc) = 

Q 2 {k) = cos(fc/2)£fcIg ] 

Q 3 (fc) = sin(fc)lgl 

Q4(fc) =sin(fc/2) gl 

Q5{k)=ek g 

Qeik) = cos(fc/2) [cr^e*^'^ ] g 

Qrik) = sin(fc) g 

Qs{k) = sin(fc/2)efc cr^g, 


(B5) 


(B6) 
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At 


FIG. 14. Evolution of the entanglement entropies for a subsystem with a number of sites £ = 1, 2, 4, 6 for the quench displayed 
in the panel and various “small” values of A. There is a perfect scaling behavior in terms of At. 


where Sk = y cos^ | + 7^ sin^ | and tan 9k = j tan |. In order to keep the corrections under control, let us isolate 
the mean-field part 


V{k, t) = I4iF(fc, ht)SV (fc, t). 


(B7) 


We obtain 


id\t5V{k,t) = hV,{k,t)SV{k,t) 

a(k.t) = ^ 


(B 8 ) 


and more explicitly 


t) = — sindfe I 




g—2'iefccr'^e^ 2 


(B9) 


where e*®'' = 3+^7 sm ^ ^ rapidly oscillating term on the second line is responsible for the irrelevance of 

y cos2 §-1-72 sin2 I 

5V in the pre-relaxation limit. Roughly speaking, one must go to second order of a perturbation theory for 5V in A 
in order to find a term proportional to the time, which is however also proportional to A^, and hence negligible in 
the pre-relaxation limit (there are no issues related to the re-exponentiation of the perturbative series). On the other 
hand, the 0(A) contributions in (B 8 ) are bounded in time, thus, at fixed At, they simply give 0(A) corrections to 
the time evolution of the observable. 

Finally, we point out that in non-interacting models the second plateau reached at large times in the pre-relaxation 
limit corresponds to the infinite time limit. 
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Appendix C: Exceptional case 7=1 

The mean-field description considered in the main text and in Ref. [69] breaks down for 7 = 1 . The reason is 
very simple: the unperturbed Hamiltonian reduces to which has infinitely many more local 

conservation laws than the XY Hamiltonian with 7 7 ^ 1. In addition, the dynamics under Hq never end up with a 
GGE, being just classical time evolution. Therefore, first of all, the mean-held description (A9) must be replaced with 

(C’)mf (0 = (^MF(i)|C>|^MF(i)) 

I^MF(t)) = [i?o + I'fMFit)) (Gl) 

I^'mf(O)) = I'I'o), 

and, second, the additional local conservation laws must be taken into account. We don’t work out this problem 
analytically but show some numerical results that testify to the emergence of a scaling regime for t ^ A“^. In Figure 
14 we plot the von Neumann entanglement entropy for subsystems with length £ = 1,2,4 ,6 after a quench from the 
Neel state to the Hamiltonian (Al) with 7 = 1, h = 0.75 and U = 0. For all subsystem sizes, the numerical data 
show a perfect scaling in terms of At at least for the interaction strength we considered, i.e. A = 0.05, 0.1, 0.2. 
Interestingly, in terms of the rescaled time and for hxed £, the different not only follow the same “average” scaling 
function, but even the amplitude of the high-frequency oscillations is modulated in such a way as to perfectly overlap: 
it seems they construct a sort of “universal” envelop. It is worth mentioning that, in this particular quench, the 
bipartite entanglement entropy (i.e. the von Neumann entropy of the semi-infinite half-chain) exhibits the same 
scaling properties. Moreover, it turns out to grow very slowly and the numerical simulations have easily reached large 
times 0 (A“^). 


